Generally, code blocks precede maps. To view code, click the “show” button.

Background and Significance

Transportation barriers pose significant challenges in providing consistent and timely care to patients. In a study using the National Health Interview Survey, nearly 21% of the United States population delayed medical care due to poor access to reliable transportation (Wolfe, McDonald, and Holmes 2020). Lack of reliable transportation can lead to missed appointments, delayed care, and inability to complete recommended treatment plans (Starbird et al. 2019; Syed, Gerber, and Sharp 2013).

Survivors of traumatic injury are a vulnerable patient population and face several barriers to accessing healthcare. Although the initial management of trauma is acute, ongoing survivorship of traumatic injuries require consistent follow up in the outpatient setting. Indeed, transportation barriers are cited in the literature as factors in influencing loss-to-follow up status (Jodoin et al. 2025). There is little insight in the literature regarding the spatial variation of transportation access and patient follow up. Understanding how transportation security varies across Chicago can give insight into how we may improve these barriers for our patient population.

Primary Means of Transportation

To understand the primary means of transportation around Chicago, data from the U.S. Census Bureau American Community Survey in 2023 was used for each census tracts aggregated over 5 years. Maps were constructed by clipping block groups into the Chicago city boundary. The census tracts and city boundary were in NAD 83 so they were reprojected to WGS 84 to make the geo-coded trauma locations and for the r5 routing engine. Block groups for the above analysis were also constructed similarly.

Population, race and ethnicity, transportation means to work, vehicle availability, and number of households below the poverty level were pulled from the 2023 ACS via the tigris package and converted to percentages.

###--------------------------means of transportation-------------------------###
#percentage of driving to work is fairly normal but the percentage who take public transit
#and the percentage who do not have a car is not; will use jenks to display
### --- plot with ggplot --- ###
#compute jenks
jenks_breaks_pdrive <- classIntervals(chi_tracts$pdrive, n = 5, style = "jenks")$brks
jenks_breaks_ppubtrans <- classIntervals(chi_tracts$ppubtrans, n = 5, style = "jenks")$brks
jenks_breaks_pwo_vehicle <- classIntervals(chi_tracts$pwo_vehicel, n = 5, style = "jenks")$brks
#categorize data
chi_tracts_jenks <- chi_tracts %>%
  mutate(drive_cat = cut(pdrive, breaks = jenks_breaks_pdrive, include.lowest = TRUE, dig.lab = 5)) %>% 
  mutate(pubtrans_cat = cut(ppubtrans, breaks = jenks_breaks_ppubtrans, include.lowest = TRUE, dig.lab = 5)) %>% 
  mutate(wo_vehicle_cat = cut(pwo_vehicel, breaks = jenks_breaks_pwo_vehicle, include.lowest = TRUE, dig.lab = 5))
#plot with jenks
#percentage who drive
drive_map <- ggplot() +
  geom_sf(
    data = chi_tracts_jenks,
    aes(fill = drive_cat),
    color = NA) +
  scale_fill_brewer(
    palette = "OrRd",
    name = "Percentage") +
  labs(
    title = "Percentage Who Drive to Work",
    caption = "Source: ACS 2023 5-year estimates (Table B08301)") +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = c(0.95, 0.80),
    legend.background = element_rect(
      fill = "white",
      color = "black"),
    legend.title = element_text(size = 12, face = "bold"),
    legend.text = element_text(size = 10),
    plot.title = element_text(
      face = "bold",
      size = 16,
      hjust = 0.5),
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank())
#percentage who do not have a car
car_map <- ggplot() +
  geom_sf(
    data = chi_tracts_jenks,
    aes(fill = wo_vehicle_cat),
    color = NA) +
  scale_fill_brewer(
    palette = "OrRd",
    name = "Percentage") +
  labs(
    title = "Percentage Who Do Not Have a Car",
    caption = "Source: ACS 2023 5-year estimates (Table B25044)") +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = c(0.95, 0.80),
    legend.background = element_rect(
      fill = "white",
      color = "black"),
    legend.title = element_text(size = 12, face = "bold"),
    legend.text = element_text(size = 10),
    plot.title = element_text(
      face = "bold",
      size = 16,
      hjust = 0.5),
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank())

#percentage who do not have a car
pubtrans_map <- ggplot() +
  geom_sf(
    data = chi_tracts_jenks,
    aes(fill = pubtrans_cat),
    color = NA) +
  scale_fill_brewer(
    palette = "OrRd",
    name = "Percentage") +
  labs(
    title = "Percentage Who Take Public Transit to Work",
    caption = "Source: ACS 2023 5-year estimates (Table B08301)") +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = c(0.95, 0.80),
    legend.background = element_rect(
      fill = "white",
      color = "black"),
    legend.title = element_text(size = 12, face = "bold"),
    legend.text = element_text(size = 10),
    plot.title = element_text(
      face = "bold",
      size = 16,
      hjust = 0.5),
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank())

#interactive versions
drive_map_int <- tm_shape(chi_tracts_jenks) +
  tm_polygons(
    col = "drive_cat",
    palette = "OrRd",
    title = "Percentage",
    border.col = NA,
    lwd = 0,
    alpha = 0.8) +
  tm_layout(
    title = "Percentage Who Drive to Work",
    bg.color = "white",
    legend.position = c("right", "top")
  )
car_map_int <- tm_shape(chi_tracts_jenks) +
  tm_polygons(
    col = "wo_vehicle_cat",
    palette = "OrRd",
    title = "Percentage",
    border.col = NA,
    lwd = 0,
    alpha = 0.8) +
  tm_layout(
    title = "Percentage Who Do Not Have a Car",
    bg.color = "white",
    legend.position = c("right", "top")
  )

pubtrans_map_int <- tm_shape(chi_tracts_jenks) +
  tm_polygons(
    col = "pubtrans_cat",
    palette = "OrRd",
    title = "Percentage",
    border.col = NA,
    lwd = 0,
    alpha = 0.8) +
  tm_layout(
    title = "Percentage Who Take Public Transit to Work",
    bg.color = "white",
    legend.position = c("right", "top")
  )

Transportation Maps

Car Availability

car_map

Drives to Work

drive_map

Public Transit to Work

pubtrans_map

Travel Time Analysis

Much of my travel time analysis was derived from the transportation chapter in Geocomputation with R (Lovelace, Nowosad, and Muenchow 2019).

Origins and Destinations

The unit of analysis chosen for travel time analysis is the census block group. Block groups have higher spatial resolution than census tracts, have fairly uniform populations (as compared to census tracts), but are more manageable for computational efforts than the highest spatial resolution - census blocks. The population centroid functions as the origin point for each block group. Population centroids approximate the geographic points of ’balance” in a block groups population density. Block groups with a population of 0 are excluded.

The destinations are level on trauma centers in Chicago as determined by the Illinois Department of Public Health. Level one trauma centers accept the highest acuity of trauma patients, and therefore will serve the population of interest - trauma patients with such sufficient injury that they need regular follow up with a trauma surgeon. I only included 9 trauma centers within or in close proximity of the city boundary by filtering the full, geo-coded trauma center list used in class. This included:

  • Advocate Christ Medical Center

  • Ascension Saint Francis

  • Mount Sinai Hospital

  • Advocate Illinois Masonic Medical Center

  • Advocate Lutheran General Hospital

  • John H. Stroger, Jr. Hospital of Cook County

  • Loyola University Medical Center

  • Northwestern Memorial Hospital

  • University of Chicago

###--------------------------origins and destinations------------------------###
#origins - centroids
USbgs_popcenter20 <- block_group2020
ILbgs_popcenter20 <- USbgs_popcenter20 %>%
  filter(STATEFP == 17) %>% # Illinois == 17
  mutate(GEOID = paste0(STATEFP, COUNTYFP, TRACTCE, BLKGRPCE)) # create a single tract GEOID similar to our other tract datasets
chi_bgs_popcenter20 <- ILbgs_popcenter20 %>%
  dplyr::select(GEOID, POPULATION, LONGITUDE, LATITUDE) %>%
  filter(GEOID %in% chi_bgs$GEOID) %>% # filter Illinois block groups on Chicago block groups only
  filter(POPULATION > 0) %>% 
  st_as_sf(coords = c("LONGITUDE", "LATITUDE"), # shapefile creation
           crs = 4326) # needs to be WGS 84
tmap_mode('view')
origins_map <- tm_shape(chi_bgs) +
  tm_borders() +
  tm_shape(chi_bgs_popcenter20) +
  tm_dots(fill = "blue",
          size = 0.25)
#destinations - trauma centers
destinations_map <- tm_shape(chiCity) +
  tm_borders() +
  tm_shape(adult_level_one) +
  tm_dots(fill = "darkgreen",
          size = 0.25)

Origins and Destinations Maps

Origins: Block Group Centroids

Destinations: Level One Trauma Centers

Routing Network Construction

The routing network was constructed from the “r5r” package in R (Pereira et al. 2021; Conway, Byrd, and Eggermond 2018; Conway, Byrd, and Van Der Linden 2017; Conway and Stewart 2019). This package uses the open source, “Rapid Realistic Routing on Real-world and Reimagined networks,” or R5, routing engine. This routing engine uses real-world street-netowrk data from OpenStreetMap and public transportation data in the form of General Transit Feed Specification (GTFS) paired with advance computational methods that account for uncertainty in public transit reliability to generate realistic estimations of public transit travel times.

Roadway and footpath data for Cook County were obtained from OpenStreetMap, and GTFS feeds for public transportation (timetables, routes, and service details) were downloaded from the CTA, Metra, and PACE websites. These inputs were used in r5r to build a multimodal routing network. Trip parameters included a maximum walking time of 30 minutes, maximum total travel time of 2 hours, and a departure window of November 19, 2025 at 09:00 ± 10 minutes. Origins were defined as population centroids of each census block group, and destinations as all level I trauma centers.

For each origin–destination pair, r5r generated multiple feasible public transit itineraries. Average travel time, wait time, walking time, number of transfers, and number of modal “legs” (walk, bus, train) were calculated by aggregating across all itineraries within the time window. Drive times were estimated separately using the road network, assuming typical traffic for the specified departure time.

Because this process produced a very large set of origin–destination combinations, each block group was ultimately assigned to a single trauma center based on the minimum drive time. The minimum drive time was chosen to approximate the trauma center that would be closest to each block group if resident was transported by EMS.

The following code depicts the above. This code is complex and computations took several hours given the magnitude of origin-destination pairs. Therefore, the code is displayed below for reference but will not run. The output is saved as an RData file for direct use in subsequent sections. A sample of the output of this file is shown below.

head(final_min_time_robust)
################################################################################
# Project: Trauma Transportation Security and Follow Up
# Script: Robust Transit Time Matrix Creation
# Purpose:
#   Create a travel time matrix accross multiple days and times using previous
#   analysis
#
# Author: Brayden Seal
# Created: 12-26-2025
# Updated: X
# Notes:
#   this uses r5r for matrix calculations which needs JDK21 downloaded
#   
################################################################################
###-------------------------------Set Up-------------------------------------###
library(matrixStats)
library(SpatialAcc)
library(leaflet)
library(leafsync)
library(tidyverse)
library(tidycensus)
library(tmap)
library(RColorBrewer)
library(rmapshaper)
library(tigris)
library(units)
library(mapboxapi)
library(sf)
library(data.table)
library(ggplot2)
library(dplyr)
library(tidytransit)
library(gtfsrouter)
library(purrr)
library(tidyr)

options(java.parameters = '-Xmx12G') #ensures appropriate memory allocated for java
#must be done prior to loading r5r package
Sys.setenv(JAVA_HOME = system("/usr/libexec/java_home -v 21", intern = TRUE)) #specify location and use of JDK 21
Sys.getenv("JAVA_HOME")
system("java -version")
library(r5r)
###--------------------------------Data--------------------------------------###
# A road network data set from OpenStreetMap in .pbf format (mandatory)
# A public transport feed in GTFS.zip format
setwd("~/Desktop/Trauma Transportation/Trauma Follow Up and Transportation Security/Robust TTM Creation")
load('base_maps.RData')
load('origins-destinations.RData')
load('origins_mindrive.RData')
mb_access_token("pk.eyJ1IjoiZ2hvc3Q0NDA4IiwiYSI6ImNtajlqd2l0bTA1eGMzc3B3dzk2YWJwYWoifQ.cWMFPBXqJQfQ0cb68L4Xug", install = T)

###--------------------------Route Construction------------------------------###
data_path <- "~/Desktop/Trauma Transportation/Trauma Follow Up and Transportation Security/Robust TTM Creation/transit_files"
#list.files(data_path)
r5r_network <- build_network(data_path = data_path)

###--------------------------------Test--------------------------------------###
#read dates
# read the GTFS feed
#gtfs_CTA <- read_gtfs('~/Desktop/GIS/transit_files/CTA.zip')
#gtfs_Metra <- read_gtfs('~/Desktop/GIS/transit_files/Metra.zip')
#gtfs_PACE <- read_gtfs('~/Desktop/GIS/transit_files/PACE.zip')

# get calendar dates (from calendar.txt)
#calendar_dates <- gtfs_CTA$calendar   

origin <- data.frame(id = '474 N Lakeshore',
                     lat = as.numeric(41.89108953920192),
                     lon = as.numeric(-87.61458686048464))

destination <- data.frame(id = 'Millenium Station',
                          lat = as.numeric(41.884589597182945),
                          lon = as.numeric(-87.6247049739245))

# set departure datetime input
mode <- c("CAR")
max_walk_time <- 30 # minutes
max_trip_duration <- 120 # minutes
departure_datetime <- as.POSIXct("16-12-2025 12:00:00",
                                 format = "%d-%m-%Y %H:%M:%S")
ttm_test <- travel_time_matrix(
  r5r_network,
  origins = origin,
  destinations = sfh,
  mode = mode,
  time_window = 60,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  max_trip_duration = max_trip_duration
)

# calculate detailed itineraries
ttm_test <- detailed_itineraries(
  r5r_network,
  origins = origin,
  destinations = destination,
  mode = mode,
  time_window = 60,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
###----------------------------Input Creations-------------------------------###
#create destinations data frame
adult_level_one$id <- adult_level_one$facility
destination <- adult_level_one %>% 
  select(id, geometry)
st_crs(destination) #needs to be WGS 84
# set departure datetime input
mode <- c("WALK", 'TRANSIT')
max_walk_time <- 30 # minutes
max_trip_duration <- 120 # minutes
departure_datetime <- as.POSIXct("18-12-2025 16:00:00",
                                 format = "%d-%m-%Y %H:%M:%S")

###---------------------Public Transit Matrix Creations----------------------###
#run this code with sequential modifications to the departure time and save after
#each run to create multiple data frames
#UCM travel time matrix
ucm <- destination %>% 
  filter(id == 'University of Chicago Emergency Department and Trauma Center')
ttm_ucm_18_16 <- detailed_itineraries(
  r5r_network,
  origins = ucm_mindrive,
  destinations = ucm,
  mode = mode,
  time_window = 60,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE)

#Christ travel time matrix
christ <- destination %>% 
  filter(id == 'Advocate Christ Medical Center')
ttm_christ_18_16 <- detailed_itineraries(
  r5r_network,
  origins = christ_mindrive,
  destinations = christ,
  mode = mode,
  time_window = 60,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
#Masonic travel time matrix
masonic <- destination %>% 
  filter(id == 'Advocate Illinois Masonic Medical Center')
ttm_masonic_18_16 <- detailed_itineraries(
  r5r_network,
  origins = masonic_mindrive,
  destinations = masonic,
  mode = mode,
  time_window = 60,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
#County travel time matrix
cook <- destination %>% 
  filter(id == 'John H. Stroger, Jr. Hospital of Cook County')
ttm_cook_18_16 <- detailed_itineraries(
  r5r_network,
  origins = cook_mindrive,
  destinations = cook,
  mode = mode,
  time_window = 60,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
#NWM travel time matrix
nmh <- destination %>% 
  filter(id == 'Northwestern Memorial Hospital')
ttm_nmh_18_16 <- detailed_itineraries(
  r5r_network,
  origins = nmh_mindrive,
  destinations = nmh,
  mode = mode,
  time_window = 60,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
#St. Francis time matrix
sfh <- destination %>% 
  filter(id == 'Ascension Saint Francis - Emergency Room')
ttm_sfh_18_16 <- detailed_itineraries(
  r5r_network,
  origins = sfh_mindrive,
  destinations = sfh,
  mode = mode,
  time_window = 60,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
#Mount Sinai time matrix
msh <- destination %>% 
  filter(id == 'Mount Sinai Hospital')
ttm_msh_18_16 <- detailed_itineraries(
  r5r_network,
  origins = msh_mindrive,
  destinations = msh,
  mode = mode,
  time_window = 60,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
#Loyola travel time matrix
loyola <- destination %>% 
  filter(id == 'Loyola University Medical Center')
ttm_loyola_18_16 <- detailed_itineraries(
  r5r_network,
  origins = loyola_mindrive,
  destinations = loyola,
  mode = mode,
  time_window = 60,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
#Lutheran travel time matrix
lutheran <- destination %>% 
  filter(id == 'Advocate Lutheran General Hospital')
ttm_lutheran_18_16 <- detailed_itineraries(
  r5r_network,
  origins = lutheran_mindrive,
  destinations = lutheran,
  mode = mode,
  time_window = 60,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)

ttm_18_16 <- bind_rows(ttm_lutheran_18_16,
          ttm_loyola_18_16,
          ttm_msh_18_16,
          ttm_sfh_18_16,
          ttm_nmh_18_16,
          ttm_cook_18_16,
          ttm_masonic_18_16,
          ttm_christ_18_16,
          ttm_ucm_18_16)
save(ttm_18_16, file = "18_16_ttm.RData")
ttm_final_raw <- bind_rows(ttm_16_8,
                           ttm_16_12,
                           ttm_16_16,
                           ttm_17_08,
                           ttm_17_12,
                           ttm_17_16,
                           ttm_18_08,
                           ttm_18_12,
                           ttm_18_16)
save(ttm_final_raw, file = "ttm_final_raw.RData")
###-----------------------Drive Time Matrix Creations------------------------###
readRenviron("~/.Renviron") #mapbox api use
mode <- c("CAR")
max_walk_time <- 30 # minutes
max_trip_duration <- 120 # minutes
departure_datetime <- as.POSIXct("18-12-2025 16:00:00",
                                 format = "%d-%m-%Y %H:%M:%S")
lutheran_mindrive_dtm <- lutheran_mindrive %>% 
  filter(id != "170311505022") #for some reason mapbox cannot calculate drive time for this block group
# 1) Build destinations + their matching origins in ONE place (paired by name)
pairs <- tibble(
  site = c("ucm","christ","masonic","cook","nmh","sfh","msh","loyola","lutheran"),
  origins = list(
    ucm_mindrive,
    christ_mindrive,
    masonic_mindrive,
    cook_mindrive,
    nmh_mindrive,
    sfh_mindrive,
    msh_mindrive,
    loyola_mindrive,
    lutheran_mindrive_dtm
  ),
  destinations = list(
    destination %>% filter(id == "University of Chicago Emergency Department and Trauma Center"),
    destination %>% filter(id == "Advocate Christ Medical Center"),
    destination %>% filter(id == "Advocate Illinois Masonic Medical Center"),
    destination %>% filter(id == "John H. Stroger, Jr. Hospital of Cook County"),
    destination %>% filter(id == "Northwestern Memorial Hospital"),
    destination %>% filter(id == "Ascension Saint Francis - Emergency Room"),
    destination %>% filter(id == "Mount Sinai Hospital"),
    destination %>% filter(id == "Loyola University Medical Center"),
    destination %>% filter(id == "Advocate Lutheran General Hospital")
  )
)

#functions for crating r5r and mapbox (traffic) drivetime matrix
r5r_drive <- function(origins_df, destinations_df) {
  detailed_itineraries(
    r5r_network,
    origins = origins_df,
    destinations = destinations_df,
    mode = mode,
    time_window = time_window,
    departure_datetime = departure_datetime,
    max_walk_time = max_walk_time,
    shortest_path = FALSE
  )
}
mapbox_drive <- function(origins_df, destinations_df) {
  mb_matrix(origins = origins_df, 
            destinations = destinations_df,
            profile = "driving-traffic",
            output = "duration")
}
#runs through each origin-destination pairing to create drive time and traffic times
dtm_long <- pairs %>%
  mutate(
    dtm = pmap(list(origins, destinations), r5r_drive),
    dtm_traffic = pmap(list(origins, destinations), mapbox_drive),
    obj_name = paste0("dtm_", site, "_")
  )
dtm_long <- dtm_long %>%
  mutate(
    dtm = map2(
      dtm,
      dtm_traffic,
      ~ {
        stopifnot(nrow(.x) == length(.y))  #ensures correct pairings
        .x %>% mutate(drive_time_traffic = as.numeric(.y))
      }))
dtm_final_raw <- bind_rows(dtm_long$dtm)
#reintroduce the block group that mapbox could not calculate
final_row <- lutheran_mindrive %>% 
  filter(id == "170311505022")
final_entry <- detailed_itineraries(
  r5r_network,
  origins = final_row,
  destinations = lutheran,
  mode = mode,
  time_window = time_window,
  departure_datetime = departure_datetime,
  max_walk_time = max_walk_time,
  shortest_path = FALSE
)
final_entry$drive_time_traffic <- 22.3
dtm_final_raw <- bind_rows(dtm_final_raw, final_entry)
save(dtm_final_raw, file = "dtm_final_raw.RData")

###-------------------------Clean Travel Time Data---------------------------###
process_ttm <- function(df) {
  
  df2 <- df %>%
    mutate(
      walk = if_else(mode == "WALK", 1, 0),
      transfer = if_else(wait > 0, 1, 0),
      bus = if_else(mode == "BUS", 1, 0),
      cta_train = if_else(mode == "SUBWAY", 1, 0),
      metra_train = if_else(mode == "RAIL", 1, 0)
    )
  
  walk_summary <- df2 %>%
    st_drop_geometry() %>%
    filter(mode == "WALK") %>%
    group_by(from_id, option) %>%
    summarise(
      total_walk_time  = sum(segment_duration, na.rm = TRUE),
      total_walk_distance = sum(distance, na.rm = TRUE),
      .groups = "drop"
    )
  
  wait_summary <- df2 %>%
    st_drop_geometry() %>%
    group_by(from_id, option) %>%
    summarise(
      wait_total = sum(wait, na.rm = TRUE),
      total_transfers = sum(transfer, na.rm = TRUE),
      total_bus_trips = sum(bus, na.rm = TRUE),
      total_cta_train_trips = sum(cta_train, na.rm = TRUE),
      total_metra_train_trips = sum(metra_train, na.rm = TRUE),
      number_of_legs = n(),
      .groups = "drop"
    )
  
  df3 <- df2 %>%
    left_join(walk_summary, by = c("from_id", "option")) %>%
    left_join(wait_summary, by = c("from_id", "option"))
  
  final <- df3 %>%
    st_drop_geometry() %>%
    group_by(from_id, to_id) %>%
    summarise(
      number_of_options = max(option),
      avg_tt = mean(total_duration, na.rm = TRUE),
      avg_td = mean(total_distance, na.rm = TRUE),
      avg_walk_time = mean(total_walk_time, na.rm = TRUE),
      avg_walk_dist = mean(total_walk_distance, na.rm = TRUE),
      avg_wait_time = mean(wait_total, na.rm = TRUE),
      avg_transfers = mean(total_transfers, na.rm = TRUE),
      avg_bus_trips = mean(total_bus_trips, na.rm = TRUE),
      avg_cta_train_trips = mean(total_cta_train_trips, na.rm = TRUE),
      avg_metra_train = mean(total_metra_train_trips, na.rm = TRUE),
      avg_legs = mean(number_of_legs, na.rm = TRUE),
      .groups = "drop"
    )
  final <- final %>% 
    rename(GEOID = from_id)
  
  return(final)
}

ttm_16_12_clean <- process_ttm(ttm_16_12)
ttm_16_16_clean <- process_ttm(ttm_16_16)
ttm_16_08_clean <- process_ttm(ttm_16_8)
ttm_17_08_clean <- process_ttm(ttm_17_08)
ttm_17_12_clean <- process_ttm(ttm_17_12)
ttm_17_16_clean <- process_ttm(ttm_17_16)
ttm_18_08_clean <- process_ttm(ttm_18_08)
ttm_18_12_clean <- process_ttm(ttm_18_12)
ttm_18_16_clean <- process_ttm(ttm_18_16)

ttm_sum <- bind_rows(ttm_16_12_clean,
          ttm_16_16_clean,
          ttm_16_08_clean,
          ttm_17_08_clean,
          ttm_17_12_clean,
          ttm_17_16_clean,
          ttm_18_08_clean,
          ttm_18_12_clean,
          ttm_18_16_clean)

ttm_clean <- ttm_sum %>% 
  group_by(GEOID, to_id) %>% 
  summarise(
    avg_number_of_options = mean(number_of_options, na.rm = TRUE),
    avg_tt = mean(avg_tt, na.rm = TRUE),
    avg_td = mean(avg_td, na.rm = TRUE),
    avg_walk_time = mean(avg_walk_time, na.rm = TRUE),
    avg_walk_dist = mean(avg_walk_dist, na.rm = TRUE),
    avg_wait_time = mean(avg_wait_time, na.rm = TRUE),
    avg_transfers = mean(avg_transfers, na.rm = TRUE),
    avg_bus_trips = mean(avg_bus_trips, na.rm = TRUE),
    avg_cta_train_trips = mean(avg_cta_train_trips, na.rm = TRUE),
    avg_metra_train = mean(avg_metra_train, na.rm = TRUE),
    avg_legs = mean(avg_legs, na.rm = TRUE))
save(ttm_clean, file = "ttm_clean.RData")

#drive time cleaning
dtm_clean <- dtm_final_raw %>%
  st_drop_geometry() %>% 
  rename(
    drive_time = total_duration,
    drive_distance = total_distance,
    GEOID = from_id) %>% 
  select(GEOID, drive_time, drive_time_traffic, drive_distance)
save(dtm_clean, file = "dtm_clean.RData")

###-------------------------Create Final Data Frame--------------------------###

final_min_time_robust <- ttm_clean %>% 
  full_join(dtm_clean, by = "GEOID")
final_min_time_robust <- chi_bgs %>% 
  left_join(final_min_time_robust, by = 'GEOID')
save(final_min_time_robust, file = "final_min_time_robust.RData")

Travel Time Chloropleths

Block group level data were aggregated to the census tract level for more intuitive presentation. Travel time discrepancy between public transit and driving was calculated as the difference between the estimated transit time and drive time. Therefore, this could also be thought of as the “time saved by driving” rather than taking public transportation. A “weighted” minimum travel time was then derived by combining transit and drive times according to household car ownership: transit time was multiplied by the proportion of households without a vehicle, and drive time by the proportion with a vehicle. This approach assumes that residents with access to a car would use it to reach the nearest trauma center, while those without a car would rely on public transit. This may better approximate the average minimum travel time to the nearest trauma center among residents in a census tract.

Chloropleths are constructed in ggplot with jenks calculated as previously described

###----------------------------Time variables--------------------------------###
final_min_time_robust <- final_min_time_robust %>%
  mutate(
    travel_disc = avg_tt - drive_time,
    weighted_travel = avg_tt * (pwo_vehicel / 100) +
      drive_time * (1 - pwo_vehicel / 100))
###-----------------Aggregate Travel Data to Tract Level---------------------###
cols_to_scale <- c('avg_tt', 'avg_td', 'avg_walk_time', 'avg_walk_dist',
                   'avg_wait_time', 'avg_transfers', 'avg_bus_trips', 'avg_cta_train_trips',
                   'avg_metra_train', 'avg_legs', 'drive_time', 'drive_distance',
                   'travel_disc', 'weighted_travel', 'drive_time_traffic')
#scale up each average by block group population
tracts_travel <- final_min_time_robust %>%
  mutate(across(all_of(cols_to_scale), ~ .x * total_popE)) 
#drop last charecter on GEOID to convert to tract GEOID
tracts_travel$GEOID <- substr(tracts_travel$GEOID, 1, nchar(tracts_travel$GEOID) - 1)
#summate scaled data
tracts_travel <- tracts_travel %>%
  st_drop_geometry() %>%
  group_by(GEOID) %>%
  summarise(
    total_popE = sum(total_popE, na.rm = TRUE),
    avg_tt = sum(avg_tt, na.rm = TRUE),
    avg_td = sum(avg_td, na.rm = TRUE),
    avg_walk_time = sum(avg_walk_time, na.rm = TRUE),
    avg_walk_dist = sum(avg_walk_dist, na.rm = TRUE),
    avg_wait_time = sum(avg_wait_time, na.rm = TRUE),
    avg_transfers = sum(avg_transfers, na.rm = TRUE),
    avg_bus_trips = sum(avg_bus_trips, na.rm = TRUE),
    avg_cta_train_trips = sum(avg_cta_train_trips, na.rm = TRUE),
    avg_metra_train = sum(avg_metra_train, na.rm = TRUE),
    avg_legs = sum(avg_legs, na.rm = TRUE),
    drive_time = sum(drive_time, na.rm = TRUE),
    drive_time_traffic = sum(drive_time_traffic, na.rm = TRUE),
    drive_distance = sum(drive_distance, na.rm = TRUE),
    travel_disc = sum(travel_disc, na.rm = TRUE),
    weighted_travel = sum(weighted_travel, na.rm = TRUE),
    .groups = "drop")
#scale back down to average by dividing by total population
tracts_travel <- tracts_travel %>%
  mutate(across(all_of(cols_to_scale), ~ .x / total_popE)) 
#combine with shapefile
final_tracts_robust <- chi_tracts %>% 
  left_join(tracts_travel, by = 'GEOID')
###-----------------------Travel Time Chloropleths---------------------------###
#compute jenks
jenks_breaks_transit <- classIntervals(final_tracts_robust$avg_tt, n = 5, style = "jenks")$brks
jenks_breaks_drive <- classIntervals(final_tracts_robust$drive_time, n = 5, style = "jenks")$brks
jenks_breaks_drive_traffic <- classIntervals(final_tracts_robust$drive_time_traffic, n = 5, style = "jenks")$brks
jenks_breaks_disc <- classIntervals(final_tracts_robust$travel_disc, n = 5, style = "jenks")$brks
jenks_breaks_weight <- classIntervals(final_tracts_robust$weighted_travel, n = 5, style = "jenks")$brks
#categorize data
travel_jenks <- final_tracts_robust %>%
  mutate(transittime_cat = cut(avg_tt, breaks = jenks_breaks_transit, include.lowest = TRUE, dig.lab = 5)) %>% 
  mutate(drivetime_cat = cut(drive_time, breaks = jenks_breaks_drive, include.lowest = TRUE, dig.lab = 5)) %>% 
  mutate(drivetime_traffic_cat = cut(drive_time_traffic, breaks = jenks_breaks_drive_traffic, include.lowest = TRUE, dig.lab = 5)) %>% 
  mutate(disc_cat = cut(travel_disc, breaks = jenks_breaks_disc, include.lowest = TRUE, dig.lab = 5)) %>% 
  mutate(weight_cat = cut(weighted_travel, breaks = jenks_breaks_weight, include.lowest = TRUE, dig.lab = 5)) 
#Cleaned labels for legend
pretty_interval_labels <- function(x) {
  out <- vapply(x, function(s) {
    # keep NA/blank labels as-is (so ggplot can show "NA" normally)
    if (is.na(s) || s == "" || s == "NA") return("N/A")
    
    nums <- str_extract_all(s, "-?\\d+\\.?\\d*")[[1]]
    if (length(nums) >= 2) {
      paste0(round(as.numeric(nums[1])), "–", round(as.numeric(nums[2])))
    } else {
      s 
    }
  }, character(1))
  
  out}

Travel Time Maps

Public Transit

#transit time
transittime_map <- ggplot() +
  geom_sf(
    data = travel_jenks,
    aes(fill = transittime_cat),
    color = NA) +
  scale_fill_brewer(
    palette = "OrRd",
    name = "Travel Time (min)",
    labels = pretty_interval_labels) +
  geom_sf(
    data = adult_level_one,
    color = "darkred",
    fill = "darkred",
    size = 2,
    shape = 21,
    stroke = 0.8) +
  labs(
    title = "Estimated Public Transit Travel Times") +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = c(0.05, 0.25),
    legend.background = element_rect(
      fill = "white",
      color = "black"),
    legend.title = element_text(size = 12, face = "bold"),
    legend.text = element_text(size = 10),
    plot.title = element_text(
      face = "bold",
      size = 16,
      hjust = 0.5),
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank())
transittime_map

Driving

#drive time
drivetime_map <- ggplot() +
  geom_sf(
    data = travel_jenks,
    aes(fill = drivetime_cat),
    color = NA) +
  scale_fill_brewer(
    palette = "OrRd",
    name = "Travel Time (min)",
    labels = pretty_interval_labels) +
  geom_sf(
    data = adult_level_one,
    color = "darkred",
    fill = "darkred",
    size = 2,
    shape = 21,
    stroke = 0.8) +
  labs(
    title = "Estimated Driving Travel Times") +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = c(0.05, 0.25),
    legend.background = element_rect(
      fill = "white",
      color = "black"),
    legend.title = element_text(size = 12, face = "bold"),
    legend.text = element_text(size = 10),
    plot.title = element_text(
      face = "bold",
      size = 16,
      hjust = 0.5),
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank())
drivetime_map

Travel Discrepancy

#travel time discrepency time
disc_map <- ggplot() +
  geom_sf(
    data = travel_jenks,
    aes(fill = disc_cat),
    color = NA) +
  scale_fill_brewer(
    palette = "OrRd",
    name = "Travel Time (min)",
    labels = pretty_interval_labels) +
  geom_sf(
    data = adult_level_one,
    color = "darkred",
    fill = "darkred",
    size = 2,
    shape = 21,
    stroke = 0.8) +
  labs(
    title = "Travel Time Saved by Driving") +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = c(0.05, 0.25),
    legend.background = element_rect(
      fill = "white",
      color = "black"),
    legend.title = element_text(size = 12, face = "bold"),
    legend.text = element_text(size = 10),
    plot.title = element_text(
      face = "bold",
      size = 16,
      hjust = 0.5),
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank())
disc_map

Weighted Travel Time

#weighted travel times
weight_map <- ggplot() +
  geom_sf(
    data = travel_jenks,
    aes(fill = weight_cat),
    color = NA) +
  scale_fill_brewer(
    palette = "OrRd",
    name = "Travel Time (min)",
    labels = pretty_interval_labels) +
  geom_sf(
    data = adult_level_one,
    color = "darkred",
    fill = "darkred",
    size = 2,
    shape = 21,
    stroke = 0.8) +
  labs(
    title = "Weighted Travel Times") +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = c(0.05, 0.25),
    legend.background = element_rect(
      fill = "white",
      color = "black"),
    legend.title = element_text(size = 12, face = "bold"),
    legend.text = element_text(size = 10),
    plot.title = element_text(
      face = "bold",
      size = 16,
      hjust = 0.5),
    panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank())
weight_map

Spatial Autocorrelation

Spatial clustering of extremes in drive times, public transit times, and weighted travel times was assessed via spatial autocorrelation. Global Moran’s I determined overall spatial correlation while local Moran’s I identified clusters of high and low spatial correlation. Briefly, Moran’s I statistic is the weighted average similarity in the measurement of interest between some observation and all it’s nearest neighbors divided by the sum of the sample variance. Nearest neighbors are determined by queen’s contiguity and measures of similarity via spatial lag with row-standardized weights (i.e. the average of the measurement of interest across all neighbors for a single observation).

Three census block groups with no population were excluded from prior analysis but will need to be included for spatial autocorrelation to avoid statistical artifact. Drive times and transit times for the centroids of each block group excluded were calculated and included in spatial autocorrelation measures. Given the populations of these block groups are zero, there is no data on household car ownership and therefore weighted travel times cannot be calculated. Their weighted travel times are approximated to be the mean weighted travel time.

Global Moran statistics for drive times, transit times, and weighted travel times all indicated spatial autocorrelation (pseudo-p < 0.01). Local Moran’s I identified clustered areas of high access (travel times less than one standard deviation from the mean and a Moran’s I > 0) and areas of low access (travel times greater than one standard deviation from the mean and a Moran’s I > 0). Scatter plots and maps below. Note that there was one census tract with significant dispersion in the weighted travel times where the tract experiences low access but is surrounded by high access.

#some tracts have no drive time because they were excluded from the routing network
#as they had no population - O'Hare, Midway, and a small Railroad tract
#these tracts are replaced with manually computed travel times
spatial_final_tracts <- final_tracts_robust
#O'Hare
spatial_final_tracts[ spatial_final_tracts$GEOID == 17031980100, "drive_time"] <- 16.67
spatial_final_tracts[ spatial_final_tracts$GEOID == 17031980100, "avg_tt"] <- 54.77
spatial_final_tracts[ spatial_final_tracts$GEOID == 17031980100, "weighted_travel"] <- mean(
  spatial_final_tracts$weighted_travel, na.rm = T)
#Railway Tract
spatial_final_tracts[ spatial_final_tracts$GEOID == 17031381700, "drive_time"] <- 10
spatial_final_tracts[ spatial_final_tracts$GEOID == 17031381700, "avg_tt"] <- 30.55
spatial_final_tracts[ spatial_final_tracts$GEOID == 17031381700, "weighted_travel"] <- mean(
  spatial_final_tracts$weighted_travel, na.rm = T)
#Midway
spatial_final_tracts[ spatial_final_tracts$GEOID == 17031980000, "drive_time"] <- 14.33
spatial_final_tracts[ spatial_final_tracts$GEOID == 17031980000, "avg_tt"] <- 37
spatial_final_tracts[ spatial_final_tracts$GEOID == 17031980000, "weighted_travel"] <- mean(
  spatial_final_tracts$weighted_travel, na.rm = T)
###-------------------------Spatial Neighbors--------------------------------###
#queen contiguity
chi_tracts_q <- poly2nb(chi_tracts, queen = TRUE) #average 6.55 links
##spatial weights matrix
#row-standardized weights for spatial neighbors
tractweights_q <- nb2listw(chi_tracts_q, style = "W")

Moran’s I and Local Moran’s I Calculations

Drive Time

#spatial lag for drive time
spatial_final_tracts$lag_drive_q <- lag.listw(tractweights_q, spatial_final_tracts$drive_time)
#plot spatial lag for drive time
ggplot(spatial_final_tracts, aes(x = drive_time, y = lag_drive_q)) + 
  geom_point(alpha = 0.3) + 
  geom_abline(color = "red") + 
  theme_minimal() + 
  labs(title = "Estimated Drive Time (min) to nearest trauma center by census tract, Chicago",
       x = "Drive Time (min)",
       y = "Spatial lag, Drive Time (min)", 
       caption = "Spatial relationships based on queens-case polygon contiguity.")

### --- standardized data --- ###
mean_drivetime = mean(spatial_final_tracts$drive_time)
sd_drivetime = sqrt(var(spatial_final_tracts$drive_time))
mean_lag_drive_q = mean(spatial_final_tracts$lag_drive_q)
sd_lag_drive_q = sqrt(var(spatial_final_tracts$lag_drive_q))
spatial_final_tracts$std_drivetime <- (spatial_final_tracts$drive_time - mean_drivetime)/sd_drivetime
spatial_final_tracts$std_lag_drive_q <- (spatial_final_tracts$lag_drive_q - mean_lag_drive_q)/sd_lag_drive_q
#map of high stadard deviation census tracts
tm_shape(spatial_final_tracts) +
  tm_borders(col = "black") +
  tm_shape(spatial_final_tracts[(spatial_final_tracts$std_drivetime > 2 & spatial_final_tracts$std_lag_drive_q > 1),]) +
  tm_polygons(col = "red",
              border.col = "black") +
  tm_layout(main.title = "Drive Time High Std Deviation",
            main.title.size = 0.7,
            frame = F,
            legend.frame = T,
            legend.outside = F,
            legend.title.size = 0.7,
            legend.text.size = 0.5)
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
### --- Moran's I --- ###
set.seed(1234)
queen_rwstd_drive <- moran.mc(spatial_final_tracts$drive_time,
                             listw = tractweights_q, 
                             nsim = 999) 
queen_rwstd_drive
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  spatial_final_tracts$drive_time 
## weights: tractweights_q  
## number of simulations + 1: 1000 
## 
## statistic = 0.9185, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater
moran.plot(spatial_final_tracts$drive_time, listw=tractweights_q)

### --- Local Moran's I --- ###
set.seed(1234)
local_drive <- localmoran_perm(
  spatial_final_tracts$std_drivetime, 
  tractweights_q, 
  nsim = 999L, 
  alternative = "two.sided"
) %>%
  as_tibble() %>%
  set_names(c("local_i", "exp_i", "var_i", "z_i", "p_i",
              "p_i_sim", "pi_sim_folded", "skewness", "kurtosis"))
#data frame with local I
local_drive_df <- spatial_final_tracts %>%
  select(GEOID, std_drivetime,std_lag_drive_q) %>%
  bind_cols(local_drive)
#clusters data frame
drive_clusters <- local_drive_df %>%
  mutate(local_cluster = case_when(
    p_i >= 0.01 ~ "Not significant",
    std_drivetime > 0 & local_i > 0 ~ "Low access",
    std_drivetime > 0 & local_i < 0 ~ "Low-high",
    std_drivetime < 0 & local_i > 0 ~ "High access",
    std_drivetime < 0 & local_i < 0 ~ "High-low"))
#plot local Moran's I
color_values <- c(`Low access` = "red", 
                  `Low-high` = "pink", 
                  `High access` = "blue", 
                  `High-low` = "lightblue", 
                  `Not significant` = "white")

ggplot(drive_clusters, aes(x = std_drivetime, 
                                      y = std_lag_drive_q,
                                      fill = local_cluster)) + 
  geom_point(color = "black", shape = 21, size = 2) + 
  theme_minimal() + 
  geom_hline(yintercept = 0, linetype = "dashed") + 
  geom_vline(xintercept = 0, linetype = "dashed") + 
  scale_fill_manual(values = color_values) + 
  labs(x = "Drive time (z-score)",
       y = "Spatial lag of drive time (z-score)",
       fill = "Cluster type")

#mapping access by car
drive_access_map <- ggplot() +
  geom_sf(
    data = drive_clusters,
    aes(fill = local_cluster),
    color = NA
  ) +
  geom_sf(
    data = chi_tracts,
    fill = NA,
    color = "black",
    linewidth = 0.3
  ) +
  scale_fill_manual(
    values = color_values,
    name = "Level of Access",
    na.value = "grey90"
  ) +
  labs(
    title = "Driving Accessibility"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = c(0.01, 0.05),
    legend.justification = c(0, 0),
    legend.background = element_rect(fill = "white", color = "black"),
    legend.title = element_text(size = 12, face = "bold"),
    legend.text = element_text(size = 10),
    plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank()
  )

Transit Time

###-------------------------Transit time Analysis------------------------------###

#spatial lag for transit time
spatial_final_tracts$lag_transit_q <- lag.listw(tractweights_q, spatial_final_tracts$avg_tt)
#plot spatial lag for transit time
ggplot(spatial_final_tracts, aes(x = avg_tt, y = lag_transit_q)) + 
  geom_point(alpha = 0.3) + 
  geom_abline(color = "red") + 
  theme_minimal() + 
  labs(title = "Estimated Transit Time (min) to nearest trauma center by census tract, Chicago",
       x = "Transit Time (min)",
       y = "Spatial lag, Transit Time (min)", 
       caption = "Spatial relationships based on queens-case polygon contiguity.")

### --- standardized data --- ###
mean_transit = mean(spatial_final_tracts$avg_tt)
sd_transit = sqrt(var(spatial_final_tracts$avg_tt))
mean_lag_transit_q = mean(spatial_final_tracts$lag_transit_q)
sd_lag_transit_q = sqrt(var(spatial_final_tracts$lag_transit_q))
spatial_final_tracts$std_transit <- (spatial_final_tracts$avg_tt - mean_transit)/sd_transit
spatial_final_tracts$std_lag_transit_q <- (spatial_final_tracts$lag_transit_q - mean_lag_transit_q)/sd_lag_transit_q
#map of high stadard deviation census tracts
tm_shape(spatial_final_tracts) +
  tm_borders(col = "black") +
  tm_shape(spatial_final_tracts[(spatial_final_tracts$std_transit > 2 & spatial_final_tracts$std_lag_transit_q > 1),]) +
  tm_polygons(col = "red",
              border.col = "black") +
  tm_layout(main.title = "Transit Time High Std Deviation",
            main.title.size = 0.7,
            frame = F,
            legend.frame = T,
            legend.outside = F,
            legend.title.size = 0.7,
            legend.text.size = 0.5)
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
### --- Moran's I --- ###
set.seed(1235)
queen_rwstd_transit <- moran.mc(spatial_final_tracts$avg_tt,
                              listw = tractweights_q, 
                              nsim = 999) 
queen_rwstd_transit
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  spatial_final_tracts$avg_tt 
## weights: tractweights_q  
## number of simulations + 1: 1000 
## 
## statistic = 0.89863, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater
moran.plot(spatial_final_tracts$avg_tt, listw=tractweights_q)

### --- Local Moran's I --- ###
set.seed(1235)
local_transit <- localmoran_perm(
  spatial_final_tracts$std_transit, 
  tractweights_q, 
  nsim = 999L, 
  alternative = "two.sided"
) %>%
  as_tibble() %>%
  set_names(c("local_i", "exp_i", "var_i", "z_i", "p_i",
              "p_i_sim", "pi_sim_folded", "skewness", "kurtosis"))
#data frame with local I
local_transit_df <- spatial_final_tracts %>%
  select(GEOID, std_transit,std_lag_transit_q) %>%
  bind_cols(local_transit)
#clusters data frame
transit_clusters <- local_transit_df %>%
  mutate(local_cluster = case_when(
    p_i >= 0.01 ~ "Not significant",
    std_transit > 0 & local_i > 0 ~ "Low access",
    std_transit > 0 & local_i < 0 ~ "Low-high",
    std_transit < 0 & local_i > 0 ~ "High access",
    std_transit < 0 & local_i < 0 ~ "High-low"))
#plot local Moran's I
ggplot(transit_clusters, aes(x = std_transit, 
                           y = std_lag_transit_q,
                           fill = local_cluster)) + 
  geom_point(color = "black", shape = 21, size = 2) + 
  theme_minimal() + 
  geom_hline(yintercept = 0, linetype = "dashed") + 
  geom_vline(xintercept = 0, linetype = "dashed") + 
  scale_fill_manual(values = color_values) + 
  labs(x = "Transit Time (z-score)",
       y = "Spatial lag of transit time (z-score)",
       fill = "Cluster type")

#mapping access by public transit
transit_access_map <- ggplot() +
  geom_sf(
    data = transit_clusters,
    aes(fill = local_cluster),
    color = NA
  ) +
  geom_sf(
    data = chi_tracts,
    fill = NA,
    color = "black",
    linewidth = 0.3
  ) +
  scale_fill_manual(
    values = color_values,
    name = "Level of Access",
    na.value = "grey90"
  ) +
  labs(
    title = "Public Transit Accessibility"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = c(0.01, 0.05),
    legend.justification = c(0, 0),
    legend.background = element_rect(fill = "white", color = "black"),
    legend.title = element_text(size = 12, face = "bold"),
    legend.text = element_text(size = 10),
    plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank()
  )

Weighted Travel Time

###----------------------Weighted Travel Time Analysis-----------------------###

#spatial lag for wieghted travel
spatial_final_tracts$lag_weight_q <- lag.listw(tractweights_q, spatial_final_tracts$weighted_travel)
#plot spatial lag for weighted time
ggplot(spatial_final_tracts, aes(x = weighted_travel, y = lag_weight_q)) + 
  geom_point(alpha = 0.3) + 
  geom_abline(color = "red") + 
  theme_minimal() + 
  labs(title = "Weighted Travel Time (min) to nearest trauma center by census tract, Chicago",
       x = "Weighted Travel Time (min)",
       y = "Spatial lag, Weighted Travel Time (min)", 
       caption = "Spatial relationships based on queens-case polygon contiguity.")

### --- standardized data --- ###
mean_weight = mean(spatial_final_tracts$weighted_travel)
sd_weight = sqrt(var(spatial_final_tracts$weighted_travel))
mean_lag_weight_q = mean(spatial_final_tracts$lag_weight_q)
sd_lag_weight_q = sqrt(var(spatial_final_tracts$lag_weight_q))
spatial_final_tracts$std_weight <- (spatial_final_tracts$weighted_travel - mean_weight)/sd_weight
spatial_final_tracts$std_lag_weight_q <- (spatial_final_tracts$lag_weight_q - mean_lag_weight_q)/sd_lag_weight_q
#map of high stadard deviation census tracts
tm_shape(spatial_final_tracts) +
  tm_borders(col = "black") +
  tm_shape(spatial_final_tracts[(spatial_final_tracts$std_weight > 2 & spatial_final_tracts$std_lag_weight_q > 1),]) +
  tm_polygons(col = "red",
              border.col = "black") +
  tm_layout(main.title = "Weighted Travel Time High Std Deviation",
            main.title.size = 0.7,
            frame = F,
            legend.frame = T,
            legend.outside = F,
            legend.title.size = 0.7,
            legend.text.size = 0.5)
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
### --- Moran's I --- ###
set.seed(1236)
queen_rwstd_weight <- moran.mc(spatial_final_tracts$weighted_travel,
                                listw = tractweights_q, 
                                nsim = 999) 
queen_rwstd_weight
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  spatial_final_tracts$weighted_travel 
## weights: tractweights_q  
## number of simulations + 1: 1000 
## 
## statistic = 0.79538, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater
moran.plot(spatial_final_tracts$weighted_travel, listw=tractweights_q)

### --- Local Moran's I --- ###
set.seed(1236)
local_weight <- localmoran_perm(
  spatial_final_tracts$std_weight, 
  tractweights_q, 
  nsim = 999L, 
  alternative = "two.sided"
) %>%
  as_tibble() %>%
  set_names(c("local_i", "exp_i", "var_i", "z_i", "p_i",
              "p_i_sim", "pi_sim_folded", "skewness", "kurtosis"))
#data frame with local I
local_weight_df <- spatial_final_tracts %>%
  select(GEOID, std_weight,std_lag_weight_q) %>%
  bind_cols(local_weight)
#clusters data frame
weight_clusters <- local_weight_df %>%
  mutate(local_cluster = case_when(
    p_i >= 0.01 ~ "Not significant",
    std_weight > 0 & local_i > 0 ~ "Low access",
    std_weight > 0 & local_i < 0 ~ "Low-high",
    std_weight < 0 & local_i > 0 ~ "High access",
    std_weight < 0 & local_i < 0 ~ "High-low"))
#plot local Moran's I
ggplot(weight_clusters, aes(x = std_weight, 
                             y = std_lag_weight_q,
                             fill = local_cluster)) + 
  geom_point(color = "black", shape = 21, size = 2) + 
  theme_minimal() + 
  geom_hline(yintercept = 0, linetype = "dashed") + 
  geom_vline(xintercept = 0, linetype = "dashed") + 
  scale_fill_manual(values = color_values) + 
  labs(x = "Weighted Travel Time (z-score)",
       y = "Spatial lag of weighted travel time (z-score)",
       fill = "Cluster type")

#mapping access by car
#mapping access by weighted travel time
weight_access_map <- ggplot() +
  geom_sf(
    data = weight_clusters,
    aes(fill = local_cluster),
    color = NA
  ) +
  geom_sf(
    data = chi_tracts,
    fill = NA,
    color = "black",
    linewidth = 0.3
  ) +
  scale_fill_manual(
    values = color_values,
    name = "Level of Access",
    na.value = "grey90"
  ) +
  labs(
    title = "Weighted Accessibility"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = c(0.01, 0.05),
    legend.justification = c(0, 0),
    legend.background = element_rect(fill = "white", color = "black"),
    legend.title = element_text(size = 12, face = "bold"),
    legend.text = element_text(size = 10),
    plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank()
  )

Mapping Access

Drive Time

drive_access_map

Transit Time

transit_access_map

Weighted Travel Time

weight_access_map

Demographic and Access Characteristics

Data on census tract characteristics is pulled from the Chicago health atlas and ACS 23 5-year survey. Shooting and crime data is from CPD.

Census Tract Qualities/Demographics

Ambulatory Difficulty

Residents reporting difficulty walking or climbing stairs as documented by ACS 23 (5-year)

tmap_mode('plot')
amb_dif_map <- tm_shape(atlas_tracts) +
  tm_borders(col = "black",
             lwd = 1) +
  tm_fill(col = "amb_difficulty",
          style = "jenks",
          palette = "reds",
          n = 5,
          title = "% of Residents") +
  tm_layout(main.title = "Ambulatory Difficulty (% of Residents)",
            main.title.size = 1,
            legend.outside = T)
amb_dif_map

Independent Living

Residents reporting difficulty doing errands alone such as visiting a doctor’s office or shopping, documented by ACS 23 (5-year)

ind_dif_map <- tm_shape(atlas_tracts) +
  tm_borders(col = "black",
             lwd = 1) +
  tm_fill(col = "ind_difficulty",
          style = "jenks",
          palette = "reds",
          n = 5,
          title = "% of Residents") +
  tm_layout(main.title = "Independent Living Difficulty (% of Residents)",
            main.title.size = 1,
            legend.outside = T)
ind_dif_map

Median Household Income

From the ACS 23; 13 tracts with missing median income data

median_inc_map <- tm_shape(atlas_tracts) +
  tm_borders(col = "black",
             lwd = 1) +
  tm_fill(col = "med_inc",
          style = "jenks",
          palette = "blues",
          n = 5,
          title = "Dollars") +
  tm_layout(main.title = "Median Annual Household Income (Dollars)",
            main.title.size = 1,
            legend.outside = T)
median_inc_map

Social Vulnerability Index

The social vulnerability index (SVI) is a measure of the demographic and socioeconomic factors that adversely affect communities(Flanagan et al. 2011). The measure was released by the CDC in 2011 and is derived from 16 U.S. census variables reflecting socioeconomic status, household characteristics, racial and ethnic minority status, housing type, and transportation indexed as percentiles relative to the state or country.

The below map reflects the SVI from 2022 covering Illinois. Levels of vulnerability are defined by quartiles (e.g. “very high vulnerability” refers to an SVI ≥ 0.75).

Methods language:

Neighborhood social vulnerability was measured using the CDC Social Vulnerability Index (SVI), a census-tract–level composite index derived from the American Community Survey. Census tracts were categorized as high vulnerability if their overall SVI percentile was ≥0.75.

atlas_tracts <- atlas_tracts %>% 
  mutate(
    svi_quant = case_when(
      svi >= 0 &  svi < 25    ~ 0,
      svi >= 25 & svi < 50    ~ 1,
      svi >= 50 & svi < 75    ~ 2,
      svi >= 75               ~ 3),
    svi_quant = factor(
      svi_quant,
      levels = c(0, 1, 2, 3),
      labels = c(
        "Very Low Vulnerability",
        "Low Vulnerability",
        "High Vulnerability",
        "Very High Vulnerability")))
svi_map <- tm_shape(atlas_tracts) +
  tm_borders(col = "black",
             lwd = 1) +
  tm_fill(col = "svi_quant",
          palette = "-RdBu",
          title = "percentile") +
  tm_layout(main.title = "Social Vulnerability Index (2022)",
            main.title.size = 1,
            legend.outside = T)
svi_map

Walkability Index

Walkability index: The National Walkability Index is a ranking describing the liklihood of walking being used as a mode of transport. values range from 1 to 20 with 20 being the most “walkable.” The walkability index is calculated from variables in the EPA’s Smart Location Database (SLD) within a specific census tract: intersection density, proximity to transit stops, and diversity of land use (mix of employment types and residential units). Scores are aggregated across every census tract nationally and divided into 20 quantiles corresponding to a score ranging from 1 to 20.

There are 18 census tracts with missing data

walkability_map <- tm_shape(atlas_tracts) +
  tm_borders(col = "black",
             lwd = 1) +
  tm_fill(col = "walk_index",
          style = "jenks",
          palette = "blues",
          n = 5)+
  tm_layout(main.title = "Walkability Index",
            main.title.size = 1,
            legend.outside = T)
walkability_map

Transportation Burden

Transportation burden: a measure of transportation insecurity accounting for relative cost of transportation and time spent of transportation. Data is sourced from the Climate & Economic Justice Screening Tool (CEJST). The data is presented as percentiles relative to the country.

There are 18 census tracts with missing data

burden_map <- tm_shape(atlas_tracts) +
  tm_borders(col = "black",
             lwd = 1) +
  tm_fill(col = "trans_burden",
          style = "jenks",
          palette = "reds",
          title = "Percentile",
          n = 5)+
  tm_layout(main.title = "Transportation Burden",
            main.title.size = 1,
            legend.outside = T)
burden_map

Crimes Data

Data on crimes from 2019-2025 is provided by the Chicago police department. The below maps reflect (1) homicides and non-fatal shootings and (2) crimes (violent and non-violent) committed at CTA locations.

Crime data is highly right-skewed so it is presented as percentiles.

Homicides and Non-Fatal Shootings

atlas_tracts <- atlas_tracts %>% 
  mutate(
    shootings_percentile = percent_rank(shootings)*100
  )
map_shootings <- tm_shape(atlas_tracts) +
  tm_borders(col = "black",
             lwd = 1) +
  tm_fill(col = "shootings_percentile",
          style = "quantile",
          palette = "reds",
          n = 4,
          title = "Percentile") +
  tm_layout(main.title = "Chicago Shootings (percentiles)",
            main.title.size = 1,
            legend.outside = T)
map_shootings

CTA Crimes

atlas_tracts <- atlas_tracts %>% 
  mutate(
    cta_crimes_percentile = percent_rank(cta_crimes)*100
  )
cta_crimes_map <- tm_shape(atlas_tracts) +
  tm_borders(col = "black",
             lwd = 1) +
  tm_fill(col = "cta_crimes_percentile",
          style = "quantile",
          palette = "reds",
          n = 4,
          title = "Percentile") +
  tm_layout(main.title = "CTA Crimes (percentile)",
            main.title.size = 1,
            legend.outside = T)
cta_crimes_map 

References

Conway, Matthew Wigginton, Andrew Byrd, and Michael van Eggermond. 2018. “Accounting for Uncertainty and Variation in Accessibility Metrics for Public Transport Sketch Planning.” Journal of Transport and Land Use 11 (1). https://doi.org/10.5198/jtlu.2018.1074.
Conway, Matthew Wigginton, Andrew Byrd, and Marco Van Der Linden. 2017. “Evidence-Based Transit and Land Use Sketch Planning Using Interactive Accessibility Methods on Combined Schedule and Headway-Based Networks.” Transportation Research Record: Journal of the Transportation Research Board 2653 (1): 45–53. https://doi.org/10.3141/2653-06.
Conway, Matthew Wigginton, and Anson F. Stewart. 2019. “Getting Charlie Off the MTA: A Multiobjective Optimization Method to Account for Cost Constraints in Public Transit Accessibility Metrics.” International Journal of Geographical Information Science 33 (9): 1759–87. https://doi.org/10.1080/13658816.2019.1605075.
Flanagan, Barry E., Edward W. Gregory, Elaine J Hallisey, Janet L. Heitgerd, and Brian Lewis. 2011. “A Social Vulnerability Index for Disaster Management.” Journal of Homeland Security and Emergency Management 8 (1): 0000102202154773551792. https://doi.org/10.2202/1547-7355.1792.
Jodoin, Zachary, Daanish Sheikh, Cameron Atkinson, Loc Uyen Vo, Alvaro Moreira, Christina Brady, and Boris Zelle. 2025. “Paediatric ballistic fracture patients: who has poor follow-up and why?” International Orthopaedics 49 (6): 1451–60. https://doi.org/10.1007/s00264-025-06506-3.
Lovelace, Robin, Jakub Nowosad, and Jannes Muenchow. 2019. Chapter 13 Transportation | Geocomputation with R. https://r.geocompx.org/transport.html.
Pereira, Rafael H. M., Marcus Saraiva, Daniel Herszenhut, Carlos Kaue Vieira Braga, and Matthew Wigginton Conway. 2021. “R5r: Rapid Realistic Routing on Multimodal Transport Networks with R5 in R.” Findings, March. https://doi.org/10.32866/001c.21262.
Starbird, Laura E., Caitlin DiMaina, Chun-An Sun, and Hae-Ra Han. 2019. “A Systematic Review of Interventions to Minimize Transportation Barriers Among People with Chronic Diseases.” Journal of Community Health 44 (2): 400–411. https://doi.org/10.1007/s10900-018-0572-3.
Syed, Samina T., Ben S. Gerber, and Lisa K. Sharp. 2013. “Traveling Towards Disease: Transportation Barriers to Health Care Access.” Journal of Community Health 38 (5): 976–93. https://doi.org/10.1007/s10900-013-9681-1.
Wolfe, Mary K., Noreen C. McDonald, and G. Mark Holmes. 2020. “Transportation Barriers to Health Care in the United States: Findings from the National Health Interview Survey, 19972017.” American Journal of Public Health 110 (6): 815–22. https://doi.org/10.2105/AJPH.2020.305579.